function [idx,dist] = deformknn(X,Y,k,s)
normX = sqrt(sum(X.^2));
X = bsxfun(@rdivide,X,normX);
nx = size(X,2);
[d,ny] = size(Y);
side = sqrt(d);
[px,py] = meshgrid(1:side,1:side);
idx = zeros(k,ny);
dist = zeros(k,ny);
tdis = zeros(2*s+1,nx);
for i=1:ny
    if (mod(i,100)==0)
        i
    end
    img = reshape(Y(:,i),side,side);
    for j = -s:1:s
        if (j == 0)
            ti = Y(:,i);
        else
            tx = (side-2*j-1)/(side-1)*(px(:)-1)+j+1;
            ty = py(:);
            ii = find(tx>=1 & tx<=side);
            tim = zeros(side,side);
            tx = tx(ii);
            ty = ty(ii);
            tim(ii) = interp2(px,py,img,tx,ty);
            ti = tim(:);
        end
        normti = 1/sqrt(ti'*ti);
        ti = normti*ti;
        tdis(j+s+1,:) = ti'*X;
%         subplot(1,2,1);imagesc(img);
%         subplot(1,2,2);imagesc(reshape(ti,side,side));
%         drawnow;pause;
    end
    mdis = max(tdis,[],1);
    [dist(:,i),idx(:,i)] = maxk(mdis',k);
end
end